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It is shown that the lowest excitation energies of a quantum many-fermion system in the random 
phase approximation (RPA) can be obtained by minimizing an effective classical energy functional. 
The minimum can be found very efficiently using generalized Lanczos technique. Application of the 
new technique to molecular spectra allows to compute excited states at the expense comparable to 
the ground-state calculations. As an example, the first-principle RPA excitation spectrum of C60 
molecule is computed taking into account all 240 valence electrons in the full valence space of the 
molecule. The results match linear absorption experiment within percents. 



Random phase approximation is central to the theory 
of electronic excitations in molecules and in extended 
systems [[[]. More generally, it resides in the core of 
the theory of linear response of correlated many-particle 
systems, and is widely used to describe correlation ef- 
fects in the excitations of nuclei J^||, molecules ||,|]], 
semiconductor quantum wells || , quantum dots M , and 
bulk materials |p| . Large systems of RPA- type equations 
are of particular significance in photochemistry of large 
molecules. Understanding such processes as photosyn- 
thesis and light reception in vision requires detailed de- 
scription of the evolution of large biological molecular 
complexes upon optical excitation , which is governed 
by the configuration of excited-state adiabatic surfaces 
pfjf . Yet, unlike ground-state molecular calculations, 
which are now considered a routine job with the pow- 
erful program packages at hand , the modeling of the 
excited states is a much more difficult task p2| . 

The reason behind this difficulty is that electronic cor- 
relations are usually much more pronounced in the ex- 
cited states. In other words, when the ground state elec- 
tronic wavefunction can be reliably approximated using 
Hartree-Fock (HF) or density functional theory (DFT), 
the electronic wavefunction in the excited state cannot 
be described as a single Slater determinant. 

RPA is one of the standard tools to treat electronic cor- 
relations in the excited sates of quantum many-particle 
systems. It belongs to a broader family of so-called 
time- dependent techniques, such as time-dependent HF 
(TDHF) and time-dependent DFT (TDDFT). These 
methods target directly the excitation energies of the sys- 
tem, by associating these energies with the frequencies of 
oscillations of the ansatz parameters when the system is 
driven out of equilibrium. The static HF or DFT ground 
state is the best Slater determinant that gives minimum 
to a certain energy functional. The TDHF and TDDFT 
describe the time evolution of the respective Slater de- 
terminant near the equilibrium. 

An assumption is made that the time-dependent wave- 
function remains a single Slater determinant at each 
moment in time. Projection of the time-dependent 
Schroedinger equation onto the set of all Slater deter- 
minants yields a set of essentially classical Hamiltonian 



equations of motion [jl|. In linear response, when the de- 
viation from equilibrium is small, the motion represents 
small oscillations, whose frequencies are associated nat- 
urally with the excitation energies of the system. In the 
small-oscillation limit TDHF is equivalent to RPA. 

RPA excitation energies are obtained as the eigenval- 
ues of a non-Hermitian matrix [B 

(-*-') V 

The N x N symmetric matrices A and B describe 
particle-particle interaction. Their matrix elements are 
simple combinations of the two-particle interaction ma- 
trix elements of the Hamiltonian in the basis of HF or- 
bitals. 

Configuration interaction singles (CIS), which is an- 
other commonly used technique to describe correlation 
effects, is in fact an approximation to RPA, and can be 
recovered by setting B = B . 

In contrast to the static HF, where the number of equa- 
tions scales linearly with the number of particles, the size 
2N of the matrix (Q) grows quadratically with the size of 
the single-particle Hilbert space. This impedes diagonal- 
ization of the matrix (^) for relatively large systems. For 
example, RPA equations for the singlet excited states of 
C60 molecule in the full valence space basis lead to the 
matrix of the size 2N = 28, 800. On the other hand, 
complete solution of the TDHF equations is not always 
necessary. In many cases, only a few low-energy excitonic 
states are of interest. 

This amounts to computing only a few extremal eigen- 
vectors of the matrix (Q) — a task similar to a standard 
problem in quantum mechanics of computing a few low- 
energy eigenstates ip of a Hermitian matrix H. The latter 
problem can be solved very efficiently using Hermitian 
Lanczos algorithm E3 ], or any of its various modifica- 
tions. In essence, the algorithm builds a Krylov subspace 
K, n of the matrix H, and then finds the best approxima- 
tion to ip in K, n by minimizing the expectation value of 
the energy (tpHip) | Q . 

There exist many variations of the Lanczos algorithm 
that allow to find eigenpairs of non-Hermitian matrices 
Q. These methods, however, loose very much of the 
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performance of the Hermitian Lanczos algorithm, be- 
cause of the lack of a minimum principle. Indeed, no 
general minimum principle exists that yields eigenvalues 
of non- Hermitian matrices p5| . 

A lot of effort has been put into developing reliable 
methods for computing selected eigenvalues of RPA-type 
matrices (0) @-§2). In m the Davidson algorithm has 
been extended to solve RPA equations as a general non- 
Hermitian eigenvalue problem. In |l7],[l8| it was modi- 
fied to preserve the special paired structure of the matrix 
(|I|). Tretiak et al. have developed a density-matrix- 
spectral-moments (DSMA) algorithm based on general- 
ized sum rules for the response theory. The symplectic 
Lanczos algorithm suggested by Mei |^] and improved 
by Benner exploits the analogy between the unitary 
transformations that preserve Hermiticity and the sym- 
plectic transformations that preserve the paired struc- 
ture of (|l|). A Newton- Raphson- type iterative procedure 
has been developed in |2^| . Finally, the oblique Lanczos 
algorithm for general non-Hermitian matrices |fl5f was 
applied to the TDHF problem in [||. 

It has been majorly overlooked that, although the 
RPA-type matrix is non-Hermitian, its block paired 
structure gives it some properties similar to the Hermi- 
tian matrices. In particular, there does exist a minimum 
principle that yields the lowest positive eigenvalue of (n|). 
It was suggested by Thouless back in 1961 and reads 
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The minimum is to be taken over all A-vectors x and y. 
The minimum always exists, since the HF stability con- 
dition keeps the numerator positive Note, that the 
denominator can be arbitrarily small, and therefore the 
expression has no maximum. 

The eigenvalue equation for the matrix (|l|) can be 
transformed into the form of Hamiltonian equations of 
motion for classical oscillations by substitution T = A+B 
and K = A - B: 



Tp = ujq, Kq — cup, 



(3) 



Here vectors q and p play the role of the conjugate canon- 
ical coordinates and momenta, while K and T are the 
matrices of stiffness and kinetic coefficients respectively. 

The lowest frequency of a harmonic Hamiltonian sys- 
tem can be obtained as a minimum of its total energy 
over all phase-space configurations {p, q} normalized by 
(pq) = l: 
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Indeed, variation of (Q) with respect to p and q yields 
Hamiltonian equations of motion (0). The minimum 



principle (Q) is equivalent to the Thouless minimum prin- 
ciple (||), where x — p + q and y = p — q. 

Two terms in the right-hand side of (4) are the kinetic 
and the potential energies at the configuration {p, q}. 
Both are positive for any p and q when the equilibrium is 
stable. The positive definiteness of K and T leads to the 
positive definiteness of the matrix in (^|), and thus, to the 
HF stability condition, making all cigcnfrcqucncies real 
& 

Minimum (4) can be easily found using the generalized 
Lanczos recursion Eal 



q i+ i = P^Tpi - ctiqi - fiiqi-i) 
Pi+i = S^iKqi - JiPi - SiPi-i), 



(5a) 
(5b) 



which generates configuration space vectors (qi,Pi) that 
span the Krylov subspace of the eigenvalue problem (|^). 
When four coefficients on, Pi, ji, and Si are chosen at 
each step i to ensure (qt+iPi) = (qt+iPi-i) = = 
(pi+i<7i_i) = 0, the vectors pi, qi form a biorthogonal 
basis, (piqj) = Sij and the matrices Kij = (qiKqj) and 
Tij = (piTpj) are symmetric tridiagonal, with the only 
nonzero matrix elements Ku — on, -Ki,j-i = Ki—i t i = 
Tu = 7i, and T^_i = Tj_ M = Si. Expanding q = ^ 
and p = "Y^diPi, we arrive at the 2n x 2n eigenvalue 
problem 



Td = ujc, Kc = cod, 



(6) 



which has the same structure as (||) . The lowest positive 
eigenvalue w m i n of (^|) give the approximation to the true 
lowest frequency oj m i n . The accuracy is found to improve 
exponentially with increasing n. 

When the lowest-frequency normal mode q^\p^ is 
found, the second-lowest normal mode q^ 2 \p^ can be 
obtained by choosing initial vectors q± and p\ orthogonal 
top'- 1 ) and q^ respectively. As follows from Eq. (Q) such 
a choice causes all vectors qi and pi to remain orthogonal 
to p^ 1 ' and cp 1 ). An oblique projection can be used to 
correct for the loss of orthogonality with respect to p^ 
and q^ that may occur at large n. Namely, the neces- 
sary amounts of q^ 1 and p^ should be subtracted from 
qi and pi respectively to ensure (qiP^) — (PiQ^- 1 ') = 0. 
Higher-frequency TDHF solutions can be found one by 
one in this way. 

In order to demonstrate the power of the new tech- 
nique, it has been applied to compute the excitation spec- 
trum of the fullerene C60 molecule. Albeit the great 
attention this molecule has received in the past sev- 
eral years [^6| , no adequate correlated calculation of the 
excited-states of C60 has been reported so far. 

As noted above, CIS in the entire particle-hole config- 
uration space can be seen as an approximation to RPA, 
where the matrix B in ([!]) is neglected. Yet, large size 
of the molecle has prevented full diagonalization of CIS 
matrix in the entire space. Calculations for C60 have 
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been reported using CIS with the CNDO/S Hamiltonian 
in the truncated space of up to 1295 out of 14400 particle- 
hole configurations Q, and using TDDFT with B-P86 
functional p^] . Also a surprisingly high lowest optically- 
allowed transition energy of 5.13 eV has been reported 
using ab-initio Rettrup-type RPA calculation ]29[ . 

The technique outlined above has allowed to solve RPA 
equations in the entire valence particle-hole configura- 
tion space of the molecule (matrix size N = 2 x 14400). 
INDO/S semiempirical parameterization of the Hamil- 
tonian |30| was used, which is essentially better than 
CNDO/S, and was shown to give especially good descrip- 
tion of the excitation spectra of 7r-conjugated molecules 
at the CIS/RPA level of theory [g,§. 




Fig. 1. C60 geometry is completely determined by the 
single and double bond lengths. 

Experimental values of 1.46 and 1.455 A |2(| were cho- 
sen for single and double bond lengths respectively, which 
completely determines the geometry of the molecule (see 
Fig. 1). INDO/S Hamiltonian matrix elements were gen- 
erated using ZINDO program |pl| . Only the singlet 
states of C60 were studied. The calculation was per- 
formed on DEC Alpha 500au workstation. Solution of 
static HF equations took about 2 min. CPU time com- 
pared to about 6 min. per each excited state. 

The present results allow for the first time for the di- 
rect comparison to the experiment. As shown in Table I, 
the energies of optically-allowed transitions obtained are 
within percents from the features observed in the linear 
absorption of C60 in solution [^2| . In Rcf . |2^] the modes 
were found to have a systematic red shift by 0.35 eV. 
No systematic shift was observed in the present study. 
Almost perfect match of all transition energies to the 
features seen in linear absorption allows to resolve the 
controversy of the assignment of the lowest optically al- 
lowed transition towards the value of 2.87 eV, which is 
opposite to the conclusion of B9], 



7 
6 

> 5 
I 4 

LU 

3 
2 
1 


Fig. 2. Excitation spectrum of C60 obtained using TDHF 
with INDO/S semiempirical Hamiltonian parameterization. 
Only Ti u states have nonzero oscillator strengths. 

Fig. 2 shows the complete excitation spectrum ob- 
tained. A total of 500 singlet excited states have been 
computed. The excitation energies were found to be de- 
generate 1, 3, 4, or 5 times in accordance with multiplic- 
ities of irreducible representations of the I/j symmetry 
group. No symmetry-induced simplification of the prob- 
lem has been used. A symmetry analysis was performed 
for each mode after it was computed and the irreducible 
representation of the symmetry group was assigned. 

Table I. C60 Experimental and theoretical electonic ex- 
citation energies and experimental oscillator strengths. Ex- 
perimental values are from linear absorption in n-hexane [^2) . 
Percent values are the deviations with respect to the experi- 
ment. 



Absorption RPA 
Experiment || INDO/S 



hid, eV 


Jose 


(full space) 
Uuj, eV 


3.04 


0.015 


2.874 (5%) 


3.30 




3.505 (6%) 


3.78 


0.37 


3.782 (0%) 


4.06 


0.10 


3.924 (3%) 


4.35 




4.287 (1%) 


4.84 




5.031 (4%) 


5.46 


2.27 


5.150 (6%) 


5.88 




5.816 (1%) 






6.008 






6.078 


6.36 




6.202 (2%) 



High symmetry of the molecule causes the majority of 
states to be optically dark. Only the states of T\ u sym- 
metry may have nonzero oscillator strengths and show up 
in linear absorption [ p6[ . It seems that the abundance of 
the singlet optically dark states below the first optically 
allowed transitions is not fully realized. The present re- 
sult may, therefore, shed some light onto a controversial 
issue of an apparent anomalously fast singlet-to-triplct 
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relaxation p6[ . 

The problem could have been simplified by taking sym- 
metry considerations into account before the RPA equa- 
tions are solved. It would be, however, opposite to the 
purpose of this letter, which is to demonstrate in the 
first place the performance of the method for a complex 
problem. In particular, specific difficulties could have 
been expected from the high level of degeneracies in the 
spectrum. No problems of that kind have been noticed. 

In conclusion, a new method is proposed for solving 
RPA-type equations with the computational effort com- 
parable to that required to solve static self-consistent 
field equations for the ground state. The method al- 
lows to compute low-energy excitonic states at the level 
of theory which may be hard or impossible to achieve 
using conventional techniques. 

As suggested in Ref. ]10| , calculation of the elec- 
tronic excitation energy at various nuclear configurations 
yields effectively the excited-state adiabatic surface of the 
molecule, provided that the ground-state adiabatic sur- 
face is known. Thus, an ability to compute the excita- 
tion energy at computational expense comparable to the 
ground-state calculation can provide a long-sought op- 
portunity to perform realistic molecular-dynamics sim- 
ulations of photochemical reactions of large biological 
molecules. 
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